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ABSTRACT 

We present simulations of collapsing 100 Mq mass cores in the context of massive 
star formation. The effect of variable initial rotational and magnetic energies on the 
formation of massive stars is studied in detail. We focus on accretion rates and on the 
question under which conditions massive Keplerian discs can form in the very early 
evolutionary stage of massive protostars. For this purpose, we perform 12 simula- 
tions with different initial conditions extending over a wide range in parameter space. 
The equations of magnetohydrodynamics (MHD) are solved under the assumption of 
ideal MHD. We find that the formation of Keplerian discs in the very early stages 
is suppressed for a mass-to- flux ratio normalised to the critical value /i below 10, in 
agreement with a series of low-mass star formation simulations. This is caused by very 
efficient magnetic braking resulting in a nearly instantaneous removal of angular mo- 
mentum from the disc. For weak magnetic fields, corresponding to /i > 10, large-scale, 
centrifugally supported discs build up with radii exceeding 100 AU. A stability analy- 
sis reveals that the discs are supported against gravitation ally induced perturbations 
by the magnetic field and tend to form single stars rather than multiple objects. We 
find protostellar accretion rates of the order of a few 10~^ Mq yr~^ which, consid- 
ering the large range covered by the initial conditions, vary only by a factor of ~ 3 
between the different simulations. We attribute this fact to two competing effects of 
magnetic fields. On the one hand, magnetic braking enhances accretion by removing 
angular momentum from the disc thus lowering the centrifugal support against grav- 
ity. On the other hand, the combined effect of magnetic pressure and magnetic tension 
counteracts gravity by exerting an outward directed force on the gas in the disc thus 
reducing the accretion onto the protostars. 

Key words: hydrodynamics - magnetic fields - methods: numerical - stars: formation 
- stars: massive. 



1 INTRODUCTION 

The question of how massive stars form is still a highly de- 



bated field of research (e.g. Zinnecker Yorke||2007 ). It is 
believed that massive star formation takes place in high mass 
molecular cloud cores with masses ranging from roughly 100 
Mq up to a few 1000 Mq. Characteristic for such cores are 
sizes of few 0.1 pc and peak densities up to 10^ cm~^ (e.g. 

Buther et al. 2007). Furthermore, from observations it is 
known that the interstellar medium as a whole is magne- 
tised (see Beck 2009 for a recent overview). Also the star 
forming cloud cores partly reveal a significant magnetisa- 
tion. The importance of the magnetic field can be estimated 
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by the mass-to-flux ratio fi normalised to the critical mass- 
to-flux ratio ( jMouschovias Spitzer^^l976J : 



Me. 



(1) 



Observed mass-to-flux ratios in high-mass star forming cores 



are typically only slightly supercritical with < 5 (Falgar- 



one et al. 1 2008) |Gir art et al.|2009||Beuther et al.|2010t indi- 
cating a significant influence of magnetic fields on the star 
formation process. In magneto-hydrodynamical simulation, 
however, also higher values of /x up to ^ 20, i.e. weaker 
magnetic fields, are found (e.g. Padoan et al. 2001[ Tilley 
k Pudritz|[2Q07| . Another common feature of star formine; 
cores are their slow rotation velocities. Observed cores have 
rotational energies normalised to the gravitational energy 
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which scatter around a mean of about 0.01 (|Goodman et al. 
1993) [Pirogov et al.||2003| ) 



In the field of low-mass star formation there is a great 
number of observations of discs and large-scale outflows 
which are the basic keystones of the widely accepted disc 
accretion scenario for low-mass star formation (see, e.g. the 
reviews by |Mac Low Klessen ] [2004) |McKee Qstriker] 



2007 ) . A similar formation scenario for massive stars is sup- 



ported by a growing number of discs and bipolar outflows 



observed around high-mass protostellar objects (see Beuther 
I&: S hepherd |2005 Cesaroni et al.|[2007 for recent reviews 



Therefore, it is worthwhile to study the influence of magnetic 
fields and rotation on the formation of discs and outflows in 
the context of massive star formation with numerical simu- 
lations. 

For low-mass star forming regions (Mcore 1 Mq), 
the influence of magnetic fields on the collapse of rotating 
cloud cores, the subsequent formation and evolution of discs 
and the launching of outfiows has received extensive atten- 
Allen et al.|2003||Banerjee k Pudr itz 2006; Pr icefcJ 



tion (e.g. 



Bat e 2007; Mellon k Li 2008; Hennebelle k Fromang 2008; 
Hennebelle k Teyssier 2008 ; Hennebelle k Ciardi 2009; Duf-. 



fin k Pudritz||2009| [Machida et al.||2011t . AU authors find 
a more or less significant infiuence of magnetic fields on the 
evolution of discs surrounding the protostars. The perhaps 
most important result of these studies is that for a mass- 
to-fiux ratio of /x < 10 the formation of Keplerian discs is 
largely suppressed. This so-called magnetic braking catas- 
trophe ( [Allen et al.|2003| [Mellon k Li|2008t turned the tra- 
ditional angular momentum problem upside down: In highly 
magnetised cores magnetic braking seems to be so eflficient 
that large-scale Keplerian discs, commonly observed around 
low-mass protostars, cannot form. Low-mass star formation 
simulations also reveal a strong impact of magnetic fields on 
the fragmentation properties of discs. In particular, strong 
magnetic fields tend to suppress disc fragmentation even in 



the presence of initial density perturbations (e.g. Hosking k 



Whitworth|2004j [Machida et al. [20051 [Hennebelle Teyssier 



2008 



Duffin k Pudritz 2009). These results are in contrast 



with the observational fact that a large fraction of low-mass 



stars are binaries (e.g. Duquennoy k Mayor[[1991 ). 

Numerical studies of the infiuence of magnetic fields 
on massive star formation, however, have received atten- 



tion only recentl y ([Banerjee k Pudritz 2007 Peters et al. 



2011 Hennebelle et al. 2011). Banerjee k Pudritz (2007) 



study the very early evolution of a protostar, its surround- 
ing disc and the outfiow performing a simulation with ex- 
tremely high spatial resolution. The protostel lar evolution 
over a timescale of some 10^ yr is examined by Peters et al. 



(2011) and Hennebelle et al. (2011), the latter authors fo- 



cussing on the effect of magnetic fields and turbulence while 
Peters et al. (2011) study the interplay of magnetic fields 



and radiation. Here, we systematically study the infiuence 
of rotation and magnetic fields on the formation and ac- 
cretion history of massive stars and examine the question 
under which conditions massive Keplerian discs can form 
in an already very early stage of protostellar evolution. We 
perform a series of collapse simulations with different initial 
rotational and magnetic energies following the protostellar 
evolution over 4000 yr. Rotational and magnetic energies 
are selected in a way to cover a large range in parameter 
space in accordance with observations and numerical simu- 



lations. Thus, we are able to detect systematic dependencies 
on the initial conditions. Furthermore, the simulations serve 
as a useful guide to select representative parameter sets for 
subsequent and more detailed studies. 

The paper is organised as follows. In Section [2l the nu- 
merical techniques and the initial conditions of the runs per- 
formed are described briefiy. The results of the simulations 
are presented in Section |3] First, the time evolution of the 
discs is presented for two representative simulations. Next, 
we analyse the velocity structure and the magnetic proper- 
ties of the gas in the midplane. Afterwards, the accretion 
histories of the protostars as well as the effects of disc frag- 
mentation are examined. In Section [H the results are dis- 
cussed in a broader context and are compared to other nu- 
merical and observational studies before we summarise our 
results in Section [S] 



2 NUMERICAL TECHNIQUES AND INITIAL 
CONDITIONS 

2.1 Numerical methods 

We present 3-dimensional, magnetohydrodynamical (MHD) 
simulations of collapsing molecular cloud cores using the 
AMR code FLASH (Fryxeh et al.|2000). We solve the equa- 
tions of ideal MHD including self-gravity. The MHD-solver 
used preserves positive states and applies well for highly su- 
personic, astrophysical problems ( Bouchut et al.|20'Q7 Waa- 



gan|2009| [Bouchut et al.|2010[[Waagan et al.|2011[ ). Using 13 
levels of refinement results in a maximum spatial resolution 
of 4.7 AU. We apply a refinement criterion that guarantees 
that the Jeans length 



Aj 




(2) 



is always resolved with at least 8 grid cells. Here Cg denotes 
the sound speed and G the gravitational constant. A reso- 
lution of 4.7 AU guarantees that we can resolve the Jeans 
length at ^ 10~^^ g cm~^ corresponding to the density limit 
where molecular gas starts to get optically thick (see also 
the Appendix). In addition, a second refinement criterion 
implemented in FLASH is applied making use of the second 
derivative of the density field. This criterion guarantees that 
protostellar discs and outfiows are resolved on the highest 
level of refinement allowing us to draw detailed conclusions 
about their properties. 

We follow the collapse of the molecular cloud core with 
increasingly higher resolution until the maximum refinement 
level is reached. As soon as the density exceeds the critical 
value of 



Pcrit = 1.78 • 10 ^^gcm ^ 



(3) 



we create a sink particle (for details see Federrath et al 



2010). Subsequently, all gas with p > pcrit within a radius of 
^sink = 12.6 AU around the sink particle is accreted onto the 
sink provided that it is gravitationally bound. The density 
cut-off at pcrit guarantees that the Jeans length is resolved 
everywhere by about 8 grid cells. The magnetic field is not 
changed during the accretion process to avoid the violation 
of the divergence- free condition. Keeping the magnetic field 
unchanged is motivated also physically by the onse t of am- 
bipolar diffusion at densities of ^ 10~^^ g cm~^ (e.g. Nakano 
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et al.|2002 ). Ambipolar diffusion would allow the gas to slip 
against the field lines, hence leaving them in the environment 
of the protostar. We briefly mention that due to numerical 
diffusion the increase of the magnetic fleld strength in the 
centre ceases at some point during the simulations therefore 
avoiding an unlimited growth of the Alfvenic velocity 

(4) 



In order to follow the thermodynamical evolution of the 
gas properly, we apply a cooling routine which includes the 
effects of molecular cooling and dust cooling as well as a 
treatment for optically thick gas (see Banerjee et al. 2006| 
for details). In this routine the cooling rates are computed 
after each hydrodynamical timestep and the state variables 
are updated accordingly for all grid cells. Using a subcy- 
cling scheme ensures that during one subcycle timestep the 
thermal energy of each cell does not change by more than 
20%. We mention that the cooling is not applied to the am- 
bient medium enclosing the molecular cloud core to prevent 
it from collapsing as well. 

As soon as an outflow is launched from the protostellar 
disc, we introduce a minimum gas density threshold within 
a radius of 67 AU around the simulation centre. Within this 
radius the gas density is kept artiflcially above a minimum 
value of 1 • 10~^^ g cm~^. Whenever a cell density falls be- 
low this threshold during one timestep, mass is added to this 
particular cell until the limit of 1 • 10~^^ g cm~^ is reached. 
Hence, we avoid the hydrodynamical timestep falling to pro- 
hibitively small values which otherwise would result from 
high Alfvenic velocities caused by low densities (see Eq. |4|. 
The effects of this artiflcial density threshold will be dis- 
cussed in detail in Section [4. II 



2.2 Initial conditions 

With the numerical methods described above we are able 
to follow the collapse of a magnetised molecular cloud core 
down to a resolution of 4.7 AU. The core has a mass of 100 
Mq and a diameter of 0.25 pc surrounded by gas with a 
density 100 times lower than the density at the core edge. 
The cubic simulation box has a length of 0.75 pc with the 
core located in its centre. Initially the core is resolved by a 
grid with a spacing of 302 AU fulfilling the applied refine- 
ment criteria. To assure pressure equilibrium at the edge 
of the core, the ambient gas has a temperature of 2000 K 
which is 100 times higher than the initially constant core 
temperature of 20 K. Calculating the Jeans mass using a 
temperature of 20 K and the average density of the core 
shows that the core contains about 56 Jeans masses and is 
therefore highly gravitationally unstable. The core itself has 
a density profile 



p(r) oc r 



(5) 



Hence, the density in the centr^ and at the edge of the 
core are initially 2.3 • 10~^^ g cm" and 4.2 • 10~^^ g cm~^. 



^ To avoid unphysically high densities in the interior of the core, 
we cut off the r"-*^- ^-profile at a radius of 0.0125 pc. Within this 
radius the density distribution follows a parabola p(r) oc (1 — 



For a mean molecular weight of /imoi — 2.3 adopted in this 
work, this corresponds to a particle density of n = 6.0 • 10^ 
and 1.1 • 10^ cm"*^, respectively. The initial core mass and 
size are comparable to that of observed massive molecular 
cloud cores (e.g. Beuther et al.|200'7 and references therein), 
where masses of 100 - 1000 Mq and sizes of a few 0.1 pc 
are found. The adopted density profile is also in good agree- 
ment with observations finding exponents around -1.5 (e.g. 



Beuther et al.|t2002b, Pirogov et al.||2003| |Pirogov||2009 ) 



In addition, the initial setup described above is comparable 
to recent studies on massive star formation in mass and/or 
density profile ( Krumholz et al.|2007 2009 Girichidis et al. 
[2Q11| iHennebelie et al.||2011| [Peters et"aL |2010a|b| |2011 ) 



which allows for quantitative comparison. 

Beside the mass distribution, which is not changed in 
this work, there are two more crucial parameters affecting 
the evolution of the cloud core, namely the magnetic field 
strength and the rotational energy. Initially, all simulations 
have a magnetic field pointing into the z-direction where Bz 
declines outwards with the cylindrical radius R as 



Bz oc R~ 



(6) 



This corresponds to a constant /3piasma = P/(B^/87r) in the 
equatorial plane where P is the thermal pressure. To guaran- 
tee VB = 0, in the initial configuration Bz is constant along 
the z-direction. We note that the chosen field configuration 
is compromise between observations and a numerically feasi- 
ble setup. However, by the time the first sink particles form, 
the field has reached a self- consistent, hourglass- like con- 
figuration as observed in massive star forming regions (e.g. 
iGirart et al.||2009| [Ah^es et al.||2011t . Therefore, we expect 
the chosen initial configuration of the magnetic field not 
to change the results significantly compared to an initially 
hourglass-shaped configuration. Initially, the cores are also 
rotating rigidly with the rotation axis parallel to the mag- 
netic field. Neither the density field nor the velocity field 
have any perturbations at the beginning of the simulation 
so that a uniform and inward directed collapse can be ex- 
pected. 

We performed 12 simulations changing the initial mag- 
netic field strength and the initial rotation frequency uj. For 
comparative purposes we also performed a simulation with 
zero magnetic field and zero rotation. The initial values of all 
simulations presented here are shown in Fig. [l] As the mag- 
netic field is changing with radius. Fig. ^ shows the initial 
ratio of magnetic energy to gravitational energy 7 and the 
corresponding normalised mass-to-fiux ratio (Eq. [T]) aver- 
aged over the entire core. The initial magnetic field strength 
at the centre of the cores ranges from 100/iG up to 1 mG. 
For comparison with 7, Fig. [l] also shows the ratio of rota- 
tional energy to gravitational energy /3rot beside the angular 
frequency. Both 7 and /3rot are calculated numerically and 
are always well below 1 as the gravitational energy is sig- 
nificantly larger than the magnetic and rotational energies. 
In addition, it can be seen from Fig. [l] that magnetic and 
rotational energies extend over more than two orders of mag- 
nitude and bracket the line where 7 equals /3rot- The initial 
values and the names of all runs are listed in Table [T] In the 
nomenclature of each run the first number give the mass-to- 
fiux ratio and the second number ^^rot multiplied by a factor 
of 100. 
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0.001 0.01 0.1 



Figure 1. Initial values of the 12 simulations performed. The 
initial rotational (/5rot) and magnetic energy (7) are normalised 
to the gravitational energy. As can be seen, the starting points 
bracket the curve where rotational and magnetic energy are equal 
(black line). The second axes show the corresponding normalised 
mass-to-flux ratio and the rotation frequency lj. 



Table 2. Formation time to of the flrst sink particle, total mass 
(the mass in brackets is the mass of the most massive sink if more 
than one is formed) of all sink particles after 4000 yr (3000 yr for 
the runs 2.6-20 and 2.6-4), the corresponding time averaged total 
accretion rate and the number of sinks created. 



Run 


to 


Msink 


Mace 


-N^sinks 




(kyrj 








26-20 


15.7 


1.26 (1.03) 


3.14 


13 


26-4 


15.2 


2.08 (1.79) 


5.19 


10 


26-0.4 


15.0 


2.93 


7.32 


1 


26-0.04 


15.0 


3.39 


8.47 


1 


10-20 


15.8 


1.28 


3.19 


1 


10-4 


15.3 


2.23 


5.57 


1 


10-0.4 


15.2 


2.98 


7.46 


1 


5.2-20 


16.7 


1.78 


4.45 


1 


5.2-4 


16.2 


2.28 


5.71 


1 


5.2-0.4 


16.1 


2.55 


6.37 


1 


2.6-20 


21.3 


1.30 


4.33 


1 


2.6-4 


20.5 


1.48 


4.93 


1 


inf-0 


15.1 


3.55 


8.86 


1 



Table 1. Performed simulations with normalised initial mass-to- 
flux-ratio /i, ratio of magnetic to gravitational energy 7, magnetic 
fleld strength in the centre B, ratio of rotational to gravitational 
energy /3rot, and the corresponding angular frequency lj. 



Run 


/i 


7 


B (/iG) 




/5rot 


UJ (10-13^-1) 


26-20 


26 


1.13 • 10-3 


132 


2 


•10-1 


7.07 


26-4 


26 


1.13 • 10-3 


132 


4 


•10-2 


3.16 


26-0.4 


26 


1.13 • 10-3 


132 


4 


• 10-3 


1.00 


26-0.04 


26 


1.13 • 10-3 


132 


4 


• 10-4 


0.316 


10-20 


10.4 


7.06 • 10-3 


330 


2 


•10-1 


7.07 


10-4 


10.4 


7.06 • 10-3 


330 


4 


•10-2 


3.16 


10-0.4 


10.4 


7.06 • 10-3 


330 


4 


• 10-3 


1.00 


5.2-20 


5.2 


2.82 • 10-2 


659 


2 


•10-1 


7.07 


5.2-4 


5.2 


2.82 • 10-2 


659 


4 


•10-2 


3.16 


5.2-0.4 


5.2 


2.82 • 10-2 


659 


4 


• 10-3 


1.00 


2.6-20 


2.6 


1.13 • 10-1 


1318 


2 


•10-1 


7.07 


2.6-4 


2.6 


1.13 • 10-1 


1318 


4 


•10-2 


3.16 


inf-0 


00 

















The considered values of rotational energies coincide 
reasonably well with values from observations, ranging from 
10~^ up to 1.4 with a mean around 0.01 (Goodma n et al.| 

Furthermore, 



1993 Pirogov et al. 2003). Furthermore, observations of 
high-mass star forming regions typically reveal mass-to-flux 
ratios /i ^ 5 and magnetic field strengths between a few 100 



/xG and a 


few mG (Lai et al. 2001] Curran & Chrysosto- 


mou|2007| 


Falgarone et al.|2008||Girart et al.|2009||Beuther 


et al.pOlO 


). Here, we explore this range with a number of 



simulations (/i ^ 5.2) but also consider initial configura- 
tions with significantly weaker magnetic fields (/i > 10). We 
recognize that turbulent motions are significant for massive 
clumps (e.g. Caselli Myers|1995 ). However, since the focus 



of this work is to investigate the formation of discs in magne- 
tised environments, we will delay the inclusion of turbulence 
to a subsequent paper. 



3 RESULTS 

In this section, we present the results of our collapse simu- 
lations where we in particular focus on the evolution of the 
discs and the accretion onto the protostars. The evolution 
and properties of the outflows launched from the protostellar 
discs will be discussed in a subsequent paper. In the follow- 
ing, we mainly limit our consideration to the phase after the 
first sink particle has been formed and only shortly describe 
the initial collapse phase. After the formation of the first 
sink particle the simulations run for further 4000 yr (except 
the runs 2.6-20 and 2.6-4 which are run for 3000 yr only) 
to study the time evolution of the outflow, the formation of 
a protostellar disc and of potential secondary sink particles 
as well as the accretion history of the sink particles. The 
computational cost of all simulations presented here sum up 
to about 1 000 000 CPU-h. 

The typical timescale for the collapse of the considered 
cloud core is its free-fall time 



Tff 



371 

32Gp 



13.9 kyr 



(7) 



where we have used the central density of p — 2.3 • 10~^^ g 
cm~^. In Table[2]we list the formation time to of the first sink 
particle for all runs performed. As can be seen, to is longer 
than Tff by a factor of 1.1 to 1.4. As run inf-0 has a pro- 
longed collapse time compared to Tff as well, we performed 
a careful numerical test showing that to can be reduced by 
a few 100 yr using a higher initial resolution. However, the 
main reason for the delay of the collapse is probably not 
a numerical issue but rather the build-up of a strong pres- 
sure gradient in the centre of the core counteracting grav- 
ity. In contrast, the unexpected fact that to for run inf-0 
is somewhat longer than for the runs 26-0.4 and 26-0.04 is 
most like a numerical issue. Turning on weak magnetic field 
in a test run with no rotation results in a slightly short- 
ened collapse time compared to inf-0. This suggest that the 
somewhat longer to in run inf-0 is an intrinsic effect of the 
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numerical scheme we will not follow up further here. How- 
ever, in each simulation subset (equal /i or /3rot) there are 
physically well motivated trends in to recognisable. As can 
be inferred from Table |2] the collapse gets slowed down with 
increasing magnetic field strength and increasing amount of 
rotational energy, both counteracting gravity. Interestingly, 
the variation of the rotational energy by two orders of mag- 
nitude at a fixed magnetic field strength changes to by no 
more than roughly 0.8 kyr, i.e. by ^ 5%. In contrast, in- 
creasing the magnetic energy by a factor of 100 for fixed 
rotational energies prolongates the collapse time by roughly 
5 kyr. However, as all cores are supercritical (/i > 1) and 
have rotational energies well below the gravitational energy 
(see Figure [T]), one cannot expect the collapse time to to be 
significantly longer than the free-fall time. In the following, 
all times refer to the time elapsed since the formation of the 
first sink particle at to- 



3.1 Global disc properties 

As a first step, we analyse the time evolution of two rep- 
resentative simulations after the formation of the first sink 
particle. For this purpose, we consider run 26-4 and run 5.2- 
4 with equal initial rotational energies but different mag- 
netic field strengths. To begin with, we compare the prop- 
erties of the gas in a well defined region of interest. This 
region is described by a cylinder around the simulation cen- 
tre with a height of 47 AU above and below the midplane. 
We mention that, although this height does not reflect the 
real disc scale height, we use it for three reasons. Firstly, 
this choice provides a numerically simple calculation of the 
desired quantities. Furthermore, it guarantees that for dif- 
ferent simulations and times identical areas are compared. 
Finally, it is used as we cannot determine a well defined disc 
scale height as even in a single run the disc height varies in 
time and radial position. For example, the frequently used 



cm ^ (see Banerjee et al. 2006 for details of the applied cool- 



approximation of the scale height (e.g. Cesaroni et al.|2007| 



H 



Cs 



(8) 



varies between a few AU and ^100 AU. This is in accor- 
dance with other estimates of the disc height made the au- 
thors. Hence we argue that our choice of 47 AU, located in 
the middle of this range, is reasonable. 

In the following we will denote this region simply as the 
disc. The quantities in the disc are averaged vertically and 
azimuthally before consideration. From the results shown 
in Fig. |2] it can be inferred that in run 26-4 an accretion 
shock occurs where both, column density and temperature 
experience a sudden increase. This is also true for run 5.2- 
4, although here the jumps in density and temperature is 
somewhat smoother. In both cases, however, the shock front 
moves outwards as time evolves. In run 26-4 the density 
profile inside the accretion shock is nearly flat with values 
around a few 10^^ cm~^ corresponding to densities of > 
10~^^ g cm"*^. For run 5.2-4 the density profile is declin- 
ing outwards and seems to decrease slightly over time. In 
both cases the mass-weighted temperature in the midplane 
increases by about one order of magnitude up to a few 100 
K. 

The temperature increase in the inner region is caused 
by the gas getting optically thick at density around 10~^^ g 



ing function). Therefore, the gas in the disc looses its ability 
to cool efficiently and cannot fast enough radiate away the 
thermal energy transfered to it by compression work. This 
results in temperatures of up to a few 100 K in the inner 
disc region consistent with observational results (see e.g. the 
review of Cesaroni et al.|20"07 and references therein) . Fur- 
ther out in the disc the gas experiences a strong temperature 
increase due to shock heating. At the accretion shock kinetic 
energy is converted quite fast into thermal energy. Assum- 
ing that gas gets decelerated by ^ 1 km s~"^ (see Fig. [2]) 
and the corresponding kinetic energy is immediately trans- 
formed into thermal energy, this would result in a temper- 
ature increase of roughly 90 K. This is in accordance with 
the observed increase of about 30 K regarding to the fact 
that a good fraction of the energy is radiated away and that 
the conversion into heat does not happen all at once. 

In Fig. [3] we show the magnetic properties of the discs. 
For the calculation of and the absolute values are 
used for averaging as both components have opposite signs 
above and below the disc and hence would cancel themselves 
out. For comparative purposes we plot Bz (initially the only 
component) and either or B^ in each panel. For the com- 
parison we only consider the radial range above 10 AU as 
the region further inside is subject to resolution effects. As 
can be seen, the different components in the midplane reach 
values up to a about 1 G. Bz is somewhat larger in run 
5.2-4 than in run 26-4 due to the five times higher, initial 
field strength in run 5.2-4. In both runs, however, B^ (bot- 
tom panel of Fig. [sj is of the order of Bz. Br is created 
by the inwards drag of the magnetic field during the col- 
lapse and later during the accretion process. In contrast to 
Br , the toroidal components of both runs differ significantly 
(see top panel of Fig. [3|. In run 26-4 is the dominant 
component in the region within the accretion shock being 
larger than Bz and Br by up to one order of magnitude. 
This is due to the fast rotation of the disc winding up the 
poloidal field (see Fig.[4]and text below). In contrast, in the 
strongly magnetised run 5.2-4, B^ is mostly smaller than 
Bz which is attributed to the lower rotation velocity in this 
case (see Fig.[4|. We also mention that all components show 
signs of an accretion shock moving outwards with time in 
accordance with the density and temperature field shown in 
Fig. [2] This is caused by the tight coupling of magnetic fields 
and matter due to ideal MHD. Furthermore, a comparison 
with the thermal pressure shows that in both runs the mag- 
netic pressure in the region within the accretion shock is 
larger than or at least equal to the thermal pressure. This 
implies that the magnetic field must play a significant role in 
the evolution of the circumstellar disc and in the accretion 
history of the protostars. 

In Fig. [4] the time evolution of the rotation (vrot), radial 
{vra.d) and Keplerian velocity (^^kep) is shown for both runs. 
For run 26-4 (left panel) a rotationally supported disc builds 
up with rotation velocities close to the Keplerian velocity. 
The maximum radius, where Vrot equals v\<,ep increases over 
time consistent with the build-up of a Keplerian disc. At the 
same time the radial velocity in the inner region nearly drops 
to zero. We emphasise that considering the behaviour within 
the central 10 AU is not conclusive as we reach the resolution 
limit at such small radii. The overall situation changes dra- 
matically when considering run 5.2-4 (right panel). Here, no 
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Figure 2. Radial profile of the column density (upper panel) and mass weighted temperature (bottom panel) for run 26-4 (left panel) 
and run 5.2-4 (right panel) at five different times after the formation of first sink particle at time to- The profiles are calculated by 
averaging azimuthally and vertically over a disc with a height of 47 AU above and below the midplane. The region below 10 AU is 
subject to resolution effects. Therefore, here and in the following plots we shaded this area to guide the readers eye. 



rotationally supported disc builds up with Vrot staying sig- 
nificantly below t'kep all the time. In fact, the absolute value 
of Vrad IS nearly as high as t'kep denoting gas which is almost 
in free-fall. Nevertheless, a flat disc-like density structure 
builds up. 

The dramatically different evolution of the gas in the 
midplane of both runs is also indicated in Fig. [5] Here 
we show the integrated column density, the magnetic field 
strength and the velocity field along a slice in the midplane. 
For run 26-4 (upper panel of Fig. [s]) a well defined Keple- 
rian disc with more or less sharp boundaries develops. The 
velocity field within the disc shows a significant rotational 
component as already seen in Fig.|4] Despite the strong frag- 
mentation occurring in the inner 200 AU, the overall disc- 
like structure is maintained. In contrast, in run 5.2-4 (lower 
panel of Fig. [s]) there is no evidence for the development 
of a Keplerian disc and the column density increases more 
or less smoothly towards the centre with nearly radial infall. 
The right panel of Fig.[5]shows the z-component of the mag- 
netic field in the midplane. As can be seen, a field strength 
of up to 1 G is reached with slightly higher values for run 
5.2-4. The close coupling of magnetic fields and matter due 
to ideal MHD is especially pronounced in the top panel of 
Fig. |5] where the strong increase in the column density at r 
^ 350 AU coincides quite well with a jump in Bz. In Fig. [6] 
the coupling of gas and magnetic field over the whole den- 
sity range is demonstrated even more clearly. Apart from 



local variations, the average of Bz and |B|, calculated in 
density bins of equal size in log-space, scales roughly as p^^^ 
or slightly weaker over more than 6 order of magnitude in 
density. This is in accordance with the scaling of a magnetic 
field in case of a spherical collapse under the condition of 
ideal MHD. 

The bubble like features seen in run 5.2-4 can be ex- 
plained as follows: during the accretion process the mag- 
netic field is dragged inwards with the gas. As we consider 
ideal MHD, the field lines cannot diffuse outwards leading 
to a constant increase in the field strength in the centre. Be- 
yond a certain point in time, the magnetic pressure is strong 
enough to effectively counteract gravity and starts to push 
material outwards. Indeed, analysing the Bz component of 
the magnetic field in the bottom right panel of Fig. [5] shows 
that the outwards moving gas is associated with a strong 
Bz. This causes Bz to decrease slightly in the very centre 
whereas at larger radii the magnetic field strength increases 
over time (see right panel of Fig.js]). 

The two simulations presented here show significant dif- 
ferences although the initial magnetic field strength is varied 
by a factor of 5 only. These systematic differences also show 
up in a more or less pronounced way when comparing all 
simulations with each other. We will discuss these differ- 
ences in detail in the next section where the complete set of 
simulations with varying initial conditions is considered. 
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Figure 3. Radial profile of the toroidal component (top) and the radial component (bottom) of the magnetic field for run 26-4 (left) 
and run 5.2-4 (right) at the same times as in Fig.|2] For comparative purposes Bz (dashed lines) is shown in both panels as well. In both 
runs Br is comparable in magnitude to Bz whereas is larger than Bz for run 26-4 and smaller than or comparable to Bz for run 
5.2-4. 




Figure 4. Radial profile of the Keplerian velocity, rotation velocity and radial velocity (negative values) for run 26-4 (left) and run 5.2-4 
(right). The velocities are averaged in the same way as in Fig. [2] and shown for same snapshots except that of t = to- For run 26-4 a 
rotationally supported structure evolves while in run 5.2-4 the rotation is clearly sub-Keplerian with radial infall close to free-fall. 



3.2 Velocity structure 

Next, we focus on the effect of the initial conditions on the 
velocity field in the midplane around the central sink par- 
ticle. We omit the time evolution as the behaviour of each 
individual run is qualitatively similar to one of the two runs 
shown before and concentrate on the situation after 4000 
yr. As already shown in Section |3.1| even small changes 



in the initial configuration of the core cause characteristi- 
cal differences. The effects of rotational and magnetic en- 
ergy on the velocity structure at the end of the simulations 
can be seen in Fig. ^ Beside the radial velocity, the ratio 
'i^rot /'^^kep is shown as well. Decreasing the initial amount of 
rotational energy for runs with fixed magnetic field strength 
(see individual panels of Fig. [7| reduces the centrifugal sup- 
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Figure 5. Column density of the disc after 2000 yr (left) and 4000 yr (middle) for run 26-4 (top) and run 5.2-4 (bottom). White dots 
mark the projected positions of the sink particles, black vectors the velocity field in the midplane. In run 26-4 the disc is well defined 
with its inner region being subject to fragmentation. In contrast, in run 5.2-4 a disc with sub-Keplerian motions forms. The bubble like 
structures in this run are caused by the high magnetic pressure in the midplane driving material outwards (see also Section [XT] ) . Right: 
z-component of the magnetic field in the midplane for run 26-4 (top) and run 5.2-4 (bottom) at the end of the simulation, i.e. after 4000 
yr. The tight correlation between field strength and matter is due to the condition of ideal MHD. Note the different spatial scales in the 
upper and lower panel. 



10-^ 



run 26-4 




run 5.2-4 













p / g cm""^ 



10-^ 



Figure 6. Average of Bz (solid lines) and of the total magnetic 
field strength (dashed lines) over density for run 26-4 (black) and 
run 5.2-4 (red) after 4000 yr. The average is calculated in density 
bins of equal size in log-space. The magnetic field is strongly cou- 
pled to the gas and scales roughly as B ~ p^/^ over more than 6 
orders of magnitude in density. 



port against gravity resulting in lower values of i^rot /'^^kep 
and consistently in higher infall velocities. Additionally, the 
differences for runs with varying firot but fixed /x seem to 
decrease when the initial field strength is increased. Fur- 



thermore, comparing runs with fixed /3rot but varying field 
strength shows that the centrifugal support, i.e. t'rot/'^^kep 
decreases with increasing field strength. This effect is pro- 
nounced strongest for high initial rotational energies. 

In principle, the simulations can be divided into two 
sets depending on the initial mass-to-flux-ratio /i. In general, 
for runs with /i > 10 centrifugally supported discs develop 
whereas for /x < 10 the formation of Keplerian discs is largely 
suppressed and only sub-Keplerian discs form at this early 
stage (see Section |4.2| for a detailed comparison with other 
numerical work) . An overview over the observed dependency 
of disc formation and fragmentation on the initial conditions 
is given in Fig.js] We note that in the runs 26-0.04 and 10-0.4 
the structure formed in the midplane strongly resembles that 
of centrifugally supported discs but with rotational velocities 
well below the Keplerian velocity. Thus, we do not denote 
them as Keplerian discs. 

In the following, we qualitatively describe the track of 
a fluid particle moving along the midplane towards the cen- 
tre. For all simulations considered in this work, the fluid 
particle first gets accelerated inwards until a radius of some 
100 AU (depending on the specific simulation considered) 
is reached. At this radius the gas experiences a decelera- 
tion, meaning its infall motion slows down. This region can 



be identified as a so-called magnetic barrier ( [Mellon Li 



2008 ) and is found in all simulations. This kind of barrier 
is different from a centrifugal barrier as here the rotation 
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Figure 7. Radial profile of the velocity field for the simulations with /i = 26 (top left), 10.4 (top right), 5.2 (bottom left) and 2.6 
(bottom right) after 4000 yr (3000 yr for /i = 2.6). For better comparison between the individual simulations the rotation velocity -Urot is 
normalised to the Keplerian velocity i;kep- The plotted quantities are again averaged azimuthally and vertically in a disc with a height of 
47 AU. For higher /3rot the ratio t'rot/'^^kep is closer to unity necessary for a Keplerian disc to form. Only in case of weak magnetic fields 
(top row) the rotation velocity reaches the Keplerian velocity while it is significantly below i;kep for strong magnetic fields. As expected, 
the absolute value of Vy-^d increases with decreasing ^'rot/'^^kep due to a lower centrifugal support against gravity. 
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Figure 8. Phase diagram of magnetic field and rotational energy 
showing the results of the simulations concerning the question of 
disc formation. For < 10 no centrifugally supported discs form 
(blue triangles) as well as for the slowly rotating cores in the 
runs 26-0.04 and 10-0.4. For weak magnetic fields Keplerian discs 
(black circles) form. Two of these discs are subject to fragmen- 
tation (red squares). The black solid line shows the curve where 
rotational and magnetic energy are equal. 



velocity is well below the Keplerian velocity which would 
be necessary to balance gravity. The reason for deceleration 
at such a magnetic barrier is twofold. Firstly, an outward 
directed thermal pressure gradient accounts for a part of 
the deceleration. This can be inferred from the significant 
increase in density and temperature at the corresponding 
radius (compare Fig. |2|. Furthermore, the magnetic field 
itself slows down the infall motion via the combination of 
an outward directed magnetic pressure gradient and the ef- 
fect of magnetic tension which can be inferred from Fig. [S] 
All components of the magnetic field experience a sudden 
increase in the region of deceleration, thus resulting in out- 
ward directed magnetic forces. Both magnetic pressure and 
magnetic tension increase in strength for smaller fi with the 
latter starting to dominate for low /i. For cases with /i < 1 
not considered here, we expect the collapse of the core per- 
pendicular to the field lines to be prevented completely by 
the Lorentz force ( Mouschovias Spitzer|1976 ). 



For radii within the magnetic barrier the velocity pro- 
files start to differ significantly from each other. For weak 
magnetic fields (top panel of Fig.|7| the gas infall speed stays 
roughly constant. Due to angular momentum conservation 
this results in an increased rotation velocity which in con- 
sequence leads to another slow down of the infall motion. 
For the highest rotational energies a centrifugal barrier is 
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encountered bringing the infall completely to halt. As seen 
in the top panel of Fig. [7| the occurrence and extension of 
the centrifugal barrier depend strongly on the initial setup. 
On the other hand, an increase in centrifugal support and 
hence a slowdown of infall gives magnetic braking more time 



to operate (Mouschovias & Paleologou 1980j). This is seen 
in the sudden drop of t'rot/'^'kep for runs with fi — 10.4 caus- 
ing a loss of centrifugal support and thus again speed-up of 
the infall motion. This is also the case for runs with low /i 
(bottom panel of Fig. |8| where magnetic braking starts to 
act efficiently directly after passing the magnetic barrier so 
that no centrifugal barrier occurs any more. The innermost 
drop of t'rad (r < 10 AU) without a corresponding increase 
in Vrot/v\<,ep is certainly caused by the limiting effect of nu- 
merical resolution as this region is only marginally resolved 
by about 3 grid cells. 

3.3 Torques 

The reason for the effective magnetic braking can be seen in 
Fig. [9] showing the edge-on view of the runs 26-4 and 5.2-4 
at two different times. As can be seen, the collapse of the 
gas has dragged in the magnetic field in the midplane pro- 
ducing a long magnetic lever arm (Alle n et al.|2003t with a 
large radial component (see also Fig.|3|. This lever arm con- 
nects the outer, slowly rotating region with the inner, fast 
rotating region, hence significantly enhancing the magnetic 
braking efficiency. As shown before, in the case of strong 
magnetic fields (/i < 10) magnetic braking is so efficient 
that large centrifugally supported discs do not form at this 
early stage. This is known as the magnetic braking catas- 
trophe and was reported previously by several authors for 
ideal MHD simulations ([Allen et al.|2003|peilon Li|2008 



Hennebe lle Froman g|2008| |Hennebelle Ciard i|2009| ). For 
a better qualitative and quantitative understanding of the 
magnetic braking effect we calculate the z-component of the 
two main torques acting on the disc, i.e. the torque exerted 
by the infalling gas 



(9) 



(10) 



Tgas = - j dW ■ [pv ■ \r X v]^) 

and the torque exerted by magnetic fields 

rmag = ^ j dV[rx{{VxB)x B)]^ . 

For the calculations we integrate over a disc with a height 
of 47 AU above and below the midplane. We omit the grav- 
itational torque as it is by two to four orders of magnitude 
smaller than Tgas and Tmag due to the nearly axisymmetric 
gravitational potential. In Fig.[lO]the time evolution of Tgas 
and Tmag is shown for the runs 26-4 and 5.2-4. The torques 
are averaged azimuthally and plotted against the radius. To 
allow for a better comparison with Tgas we plot —Tmag- A 
positive T denotes a flux of angular momentum into the disc 
while for a negative r angular momentum is removed from 
it. Hence, in both runs the magnetic field is removing angu- 
lar momentum from the disc, i.e. slowing down its rotation 
(Tmag < 0) whereas the gas exerts a positive torque on the 
disc. It can be seen that the torques increase steadily with 
time (except for run 26-4 at 3000 yr). For run 26-4 Tgas is 
nearly always larger than the magnetic torque whereas in 
run 5.2-4 Tgas is roughly balanced by —Tmag from the very 



beginning. These differences show up even more clearly in 
Fig. [11] where the torques at the end of the simulations are 
shown. For runs with a weak magnetic field (top left panel of 
Fig. 11) the torque exerted by the gas is nearly everywhere 
above the magnetic torque. Thus, there is a net flux of angu- 
lar momentum into the disc leading to the observed build-up 
of centrifugally supported discs. Only for run 26-20 —Tmag 
equals or even exceeds Tgas in the inner region although only 
on a low level. This is due to the large Keplerian disc which 
has already built up in this run having only very small infall 
velocities (see top left panel of Fig. [7|. The sharp jump of 
Tgas around r = 200 AU is caused by the accretion shock at 
the edge of the disc where Vrad drops to zero (see top left 
panel of Fig. |7|. 

Analysing the different panels in Fig. [TT] it can be seen 
that at large radii the gas torques are always larger than 
—Tmag. This implies that the magnetic field has only a small 
effect on the collapse in the outer parts. In contrast, at 
smaller radii the effect of magnetic fields gets more and 
more pronounced. With decreasing /i the magnetic torque 
approaches Tgas and in particular for the runs with /i = 5.2 
and 2.6 —Tmag is very close to Tgas- This is attributed to the 
fact that an equilibrium between Tgas and Tmag is reached 
where as much angular momentum is removed by magnetic 
braking as it is added due to the gas infall. Hence, in the case 
of strong magnetic fields the net angular momentum flux is 
roughly zero preventing a Keplerian disc from forming. 

To summarise, centrifugally supported discs at very 
early times only form in simulations with fi > 10. In these 
cases magnetic braking is not strong enough to remove an- 
gular momentum at high enough rates and hence is not 
able to prevent Keplerian disc formation (see Fig. [s] for an 
overview) . 



3.4 Accretion rates 

Closely related to the velocity structure of the matter around 
the sink particles is the accretion onto the particles them- 
selves. The general behaviour of the accretion rate with vary- 
ing initial conditions can be inferred from Table |2] where we 
list the totally accreted masses and the corresponding time 
averaged accretion rates. Additionally, in Fig. [12] the accre- 
tion rates of all runs are plotted against the initial mass-to- 
flux ratio. As listed in Table [2] in neither of the runs more 
than 4 Mq are accreted during the simulation resulting in 
time averaged accretion rates around a few 10~^ Mq yr~^. 
Interestingly, the accretion rates do not vary by more than 
a factor of about 3 between the different simulations. This is 
remarkable, considering the large range in parameter space 
covered by the initial conditions (Fig.[T]). For each set of sim- 
ulations with equal /i there is also a rough correspondence 
between the accretion rate and the infall velocity shown in 
Fig. [7| As expected, higher infall motions result in higher 
accretion rates. 

For a more detailed analysis of the accretion rates, we 
consider their time evolution in Fig. |13| We mention that for 
the runs 26-20 and 26-4 where more than one sink particle 
is formed only the accretion onto the first particle created 
is shown. A more detailed analysis of the accretion history 
in the case of fragmentation will be carried out in Section 
|3.5| As can be seen in Fig. |13| there seems to be a slight 
decrease in the accretion rates by a few 10% over time. Nu- 
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Figure 9. Edge-on view of the central region for run 26-4 (left) and run 5.2-4 (right). Superposed on the column density are the velocity 
field (black vectors) and the magnetic field lines (white lines). Also shown are the regions where the toroidal magnetic field dominates 
over the poloidal field (black lines) and the region used for calculating disc properties (dark green line). As can be seen, in run 26-4 
the greatest part of the inner region is dominated by the toroidal magnetic field whereas in run 5.2-4 only smaller parts further out are 
dominated by B^. Note the different spatial scales in the left and right panel. 




Figure 10. Gas torque (solid lines) and magnetic field torque (dashed lines) exerted on the disc for run 26-4 (left) and run 5.2-4 (right). 
The torques are averaged azimuthally and shown at the same times as in Fig.|4] As Tgas is negative, we plot — Tgas for better comparison. 
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Figure 11. Gas torque (solid lines) and magnetic field torque (dashed lines) exerted on the disc for runs with /i = 26 (top left), 10.4 (top 
right), 5.2 (bottom left) and 2.6 (bottom right) at the end of the simulations. As Tgas is negative, we plot — Tgas for better comparison. 
For high magnetic field strengths the magnetic torque roughly balances the gas torque resulting in the suppression of Keplerian disc 
formation. 
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Figure 12. Accretion rates of the different simulations plotted 
against the initial mass-to-fiux-ratio fi. Equal colours denote equal 
initial rotational energies. The accretion rates seem to converge 
with decreasing /x. 



merical studies (e.g. 'Kless en|200l" Schmeja &: Klessen|2004 ) 



Andre et al.|2000| for an 



as well as observational results (see . 

overview) indicate that there is indeed a decline in the ac- 
cretion rate by orders of magnitude, although timescales for 
this decline are typically of the order of several 10^ yr and 
are therefore much longer than in our study. Additionally, 



there occur fast variations within a factor of about 2 around 
the mean value. This fast variations are caused by moder- 
ate density perturbations developing in the midplane. Each 
time a perturbation moves through the centre, it causes the 
accretion to vary around its mean. The variations are gen- 
erally rather small (with one exception in run 26-0.4) and 
would probably be smoothed out in time by viscous effects 
in the inner disc not resolved here. However, we cannot ex- 
clude that the varying accretion rates would influence the 
protostellar evolution as proposed by stellar evolution cal- 
culations (e.g. Wuchterl Klessen|2001 Baraffe & Chabrier 

^oTol. 



The only exception where the accretion rates decrease 
significantly over time are the runs with /x = 2.6 showing a 
sharp drop in the mass accretion after roughly 2500 yr. This 
is caused by the occurrence of magnetically driven bubbles 
in the midplane as shown exemplarily in the bottom panel of 
Fig. [5] However, in the other runs the accretion rates seem to 
decrease only slightly showing no sign that the magnetically 
driven outflow, which is launched from the protostellar disc 
shortly after the formation of the protostar, can significantly 
reduce mass accretion over time. Similar results in related 
work on magnetic fields in massive star formation are found 



by [Peters et al.| ( [2011[ ) and |Hennebelle et al.| ( |2011| ) as well. 
This behaviour is due to the ongoing accretion through the 
disc which is nearly unaffected by the shut-off of gas infall 
from below or above the disc due to the outflow. This can 
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be seen in Fig. [9] where we plot the velocity structure and 
magnetic field lines in a slice along the z-axis. 

As can be inferred from Table |2] and from Fig. |12| there 
are some systematical trends in the accretion rates with 
changing initial conditions, although the overall variation 
is not very large. Increasing the overall rotational support 
against gravity, i.e. /3rot, for fixed /i results in lower accre- 
tion rates. This is in agreement with the increase in t'rot/'^^kep 
and the decrease in Vrad in the surrounding disc as shown 
in Fig. |7| As already observed for the velocity structure in 
the midplane, the differences in the accretion rates for runs 
with different /3rot but fixed /i decrease with increasing mag- 
netic field strength as can be seen in Fig. |12| This is due to 
the efficient magnetic braking removing angular momentum 
at roughly the same rate as it is transported inwards (see 
Fig. 11). As a consequence, the accretion rates are roughly 



independent of the initial amount of angular momentum. An 
even further increase in the field strength would probably re- 
sult in even lower accretion rates than the 4-5 -10"^ Mq 
yr~-^ observed for /i = 2.6 due to stronger magnetic forces 
counteracting gravity. 

The accretion rates for varying /x but fixed /Srot show 
an interesting and less clear behaviour. For high initial rota- 
tional energies, i.e. /3rot = 0.20 and 0.04, respectively, the ac- 
cretion rates first increase with decreasing fi and drop again 
for the case with /i = 2.6. For /3rot = 0.004 the accretion 



rate increases from /i = 26 to /i = 10.4 before declining 
again for lower /i. We attribute this behaviour to two com- 
peting effects of the magnetic field. On the one hand, mag- 
netic fields act to enhance accretion onto the protostar by 
magnetic braking reducing the centrifugal support against 
gravity. Hence, the effect of magnetic braking alone would 
cause increasing accretion rates with decreasing /j as it is 
indeed observed for low field strengths. The second effect 
influencing the accretion rates is the Lorentz force, i.e. the 
combination of magnetic pressure and magnetic tension, in- 
duced by strongly bent field lines (see Fig. [9]). Magnetic pres- 
sure and magnetic tension counteract gravity by exerting an 
outward force on the gas resulting in reduced accretion rates. 
The strength of this effect increases with the field strength 
thus tending to lower accretion rates with lower /i. 

The combination of both effects - magnetic braking en- 
hancing accretion and the Lorentz force counteracting ac- 
cretion - results in the observed behaviour of the accretion 
rates: By increasing the magnetic field strength at a given 
/3rot up to a certain critical value an equilibrium between 
the torques acting on the disc is reached where the removal 
of angular momentum by magnetic braking balances its in- 



wards transport (see Fig. 11). Up to this value, an increase 
in the field strength is associated with increasing accretion 
rates as observed in our simulations. Further increase in the 
field strength beyond this point cannot enhance the mag- 
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netic braking efficiency any more. In fact, now tfie increase 
(decrease, if /i is considered) ieads to deciining accretion 
rates due to tfie growing strengtfi of tfie Lorentz force coun- 
teracting gravity, in accordance witfi our findings of declin- 
ing accretion rates for strong fields with /i < 5. The exact 
value of /i, where this turnover occurs, depends on the initial 
amount of rotational energy and decreases with increasing 

/^rot- 



3.5 Disc fragmentation 

As mentioned earlier, the runs 26-20 and 26-4 show rapid 
fragmentation of the protostellar disc after the first protostar 
has formed (see Fig.js]). Due to a high amount of rotational 
energy and a weak magnetic field, magnetic braking can only 
remove a small amount of angular momentum leading to the 
formation of Keplerian discs with considerable extensions of 
a few 100 AU (see top left panel of Fig. [t] for comparison). 
As the mass load onto these discs exceeds their capability 
to transport material inwards by gravitational or viscous 



torques, the discs become unstable and fragment (e.g. Krat- 
ter et al.|2010 ). At the end of the simulations, i.e. after 4000 
yr there are 10 sink particles in run 26-4 and 13 in run 26-20. 
All other simulations show no fragmentation so far although 



some of them form a Keplerian disc (compare Fig. 12). 

The accretion histories for run 26-20 and run 26-4 are 
shown in Fig. |14| Run 26-20 exhibits a very symmetric disc 
fragmentation forming pairs of protostars at roughly the 
same time and opposite positions (as seen from the cen- 
tre). Even in their further evolution each pair develops very 
similar as can be seen from the left panel of Fig. |14| where 
the lines of each pair are nearly indistinguishable. For run 
26-4 only the evolution of the second and third sink parti- 
cle is symmetric, while at later times there is no pairwise 
formation of sink particles any more due to a nonsymmetric 
evolution of the disc. We note that in both runs only the first 
sink particle created has reached a mass above 1 Mq so far 
whereas all other particles have masses well below 0.1 Mq. 
For comparative purposes we also plot the totally accreted 
mass of all sinks formed in Fig. 14 (see also Table [2]). As can 
be seen, in both runs more than 



70 of the total mass is 
accreted onto the first sink particle. 

After about 2500 yr, i.e. after the creation of further 
sink particles, there occurs a slight but nevertheless notice- 
able decrease in the accretion onto the first sink particle 



(see the red and black lines in the top left panel of Fig. 13). 
This fragmentation-induced starvation jPeters et al.|2010b|) 
is caused by the surrounding sink particles soaking up the 
infalling material. This behaviour was also observed recently 



in related work on massive star formation (jPeters et al, 
2010a| [20TT| [Giri chidis e TaLllml] ) . 



Beside the sink masses. Fig. [14] shows the calculated 
disc masses as well. For the calculation we define the disc 
as follows: Firstly, we determine the maximum cylindrical 
radius Rmax where the gas falls below a density of 5 • 10~^^ 
g cm~^. The disc mass is then defined as the total mass 
of gas within a cylinder with a height of 47 AU above and 
below the midplane and a radius of Rmax around the centre. 
As can be seen, the disc masses lie around 1 Mq and are 
therefore somewhat below the totally accreted sink masses. 
Next, we study the stability of the discs against gravi- 



tationally induced perturbations. Disc stability is described 
by the Toomre parameter ( Toomre||1964 ) defined as 



(11) 



with the epicyclic frequency sound speed Cs, surface den- 
sity S, and gravitational constant G. Gravitational instabil- 
ity sets in when Q < 1. As magnetic fields are present in the 



discs, one can define the magnetic Toomre parameter (Kim 
k Ostriker|[200l| ) 



2 , 2\l/2 



(12) 



where va is the Alfvenic velocity taking into account all 
components of the magnetic field. In the following, we anal- 
yse the stability of the midplane in the runs 26-20, 10-20, 
5.2-20 and 26-0.4 as these simulations cover a wide range of 
initial conditions. In particular, we concentrate on the sta- 
bilising effect of the magnetic field by comparing Qm and 
Q. In Fig. [15] the face-on view of the discs after 4000 yr as 
well as the radial dependence of the Toomre parameters are 
plotted. Q and Qm are calculated numerically by the az- 
imuthally averaged values of Cs, va, and S and are shown 
for several times. We neglect the contribution of secondary 
sink particles to the column density in run 26-0.4 as their 
only effect would be to lower Q and Qm in regions which are 
already unstable and fragmenting and therefore they would 
not qualitatively change the situation. The plot shows that 
within a radius of r :^ 500 AU the Toomre parameter Q is 
almost everywhere below 1 for all runs and all times con- 
sidered. Hence, in all runs fragmentation would be expected 
to occur. As this is not the case except for run 26-20, the 
Q-parameter seems not able to properly describe the stabil- 
ity properties of the discs. Taking into account the effect of 
magnetic fields, the situation displayed in the column den- 
sity plots is described much more adequately. The value of 
Qm is higher than Q by roughly a factor of 10 in all cases 
thus pointing to much more stable configurations. Never- 
theless, for run 26-20 even Qm is below 1 with one spatially 
confined exception at 2000 yr. Hence, in this case even the 
combination of thermal and magnetic pressure cannot sta- 
bilise the disc against fragmentation as shown in the upper 
left panel of Fig. |15| Fragmentation occurs in the innermost 
part (r < 50 AU) close to the first sink particle and in a 
filamentary- like ring at r :^ 200 AU. We note that the situa- 
tion is somewhat different to run 26-4 where fragmentation 
seems to occur mainly in spiral like density perturbations 
closer to the centre (see top panel of Fig. |5]). However, the 
stability analysis of run 26-4 shows that Q and Qm are below 
1 as well. 

All other runs considered in Fig. |15| show no fragmen- 
tation so far. This agrees quite well with the findings of 
Qm — 1- However, there are some situations when Qm is 
marginally smaller than 1 with the disc still showing no signs 



of fragmentation. Kratter et al. (2010), analysing the stabil 



ity of discs under purely hydrodynamical conditions, find 
that discs can be indeed stable if Q is locally smaller than 
1. The authors attribute this to the fact that Q = 1 indi- 
cates instability of axisymmetric perturbations in infinitely 
thin discs ( Toomre|1964 ). For thick discs, as it is the case in 
our simulations, they argue that the instability criterion is 



expected to be somewhat relaxed ( Goldreich & Lynden-Bell 
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Figure 14. Accretion history of run 26-20 (left) and run 26-4 (right). The dashed-dotted hues show the total mass of all sink particles. 
The formation of sink particles in run 26-20 occurs pairwise (except the first one) and also the further evolution of each pair is nearly 
indistinguishable except at the very end so that in the beginning each line represents two particles. Also shown is the evolution of the 
disc mass with time (dashed lines) which is of the order of or somewhat below the totally accreted mass. 



1965). Hence, we expect our discs to be stable for Toomre 
values even slightly smaller than 1. 

To summarise, the comparison of Q and Qm shows that 
the stabilising effect of magnetic fields is needed to effec- 
tively prevent the discs from fragmenting. For weak mag- 
netic fields (/i = 26) and high amounts of angular momen- 
tum, however, even magnetic fields cannot stabilise the discs 
against fragmentation. Furthermore, we expect that frag- 
mentation would occur for other runs with low magnetic 
field strengths (/i ^ 10) as well if the evolution of the discs 
would be followed beyond 4000 yr. This assumption is con- 
firmed by the results that the discs show a continuing growth 
due to ongoing infall and are only marginally stable (Qm — 
!)• 



4 DISCUSSION 

4.1 Numerical caveats 

The simulations considered so far were run with a maxi- 
mum resolution of 4.7 AU. To test the resolution dependency 
of our results we performed two more simulations with the 
same initial setup as in run 26-4 but with a resolution varied 
by a factor of 4 in either direction. A detailed comparison of 
the results is shown in the appendix. The results of the res- 
olution study suggest that the spatial resolution of 4.7 AU 
used in the simulations presented in this work is sufficient to 
properly follow the dynamical evolution of the protostellar 
discs and of the accretion rates. 

The bubble like features seen in the lower panel of Fig. [5] 
also occur in other runs with strong initial magnetic fields 
and are artifacts of our assumption of ideal MHD. In par- 
ticular, for runs with /i = 2.6 the bubbles show a signifi- 
cant influence on the accretion rates (bottom right panel of 



Fig. 13) reducing the accretion after 2500 yr. However, for 
runs with /x = 5.2 the dynamical influence of such bubbles 
seems to be rather limited as we cannot detect any signifi- 
cant changes in the accretion rates (see bottom left panel of 
Fig. 13). We therefore conclude that the accretion rates are 
reliable for runs with /i down to ^ 5 and up to 2500 yr even 



for the both runs with /x = 2.6. A possible solution to avoid 
the formation of bubbles like features could be the inclusion 
of the effects of non-ideal MHD such as ambipolar diffusion 
or ohmic dissipation as done in recent work (^Puffin Pu- 



dritz|2009| [Mellon k Li|2009||Dapp k Basu|2010t . However7 
Nakano et al. (|2002") argue that ohmic dissipation starts to 
act efficiently at particle densities above 10~^^ cm~^ which 
is well above the typical densities found in our discs. Hence, 
ohmic dissipation is not expected to be capable of signifi- 
cantly reducing the magnetic flux in the centre. Ambipolar 



diffusion, however, starts to act at lower densities ( [Puffin 
Pudri'tz]|2009 ) and thus might be able to reduce the field 
strength in the centre before magnetically driven bubbles 
can occur. 

As mentioned in Section [2Tj a minimum density thresh- 
old of 1 • 10~^^ g cm~^ is applied after the outflow is 
launched. The exact amount of mass added artificially dur- 
ing the 4000 yr depends on the individual simulation but 
never exceeds a few 10~^ Mq corresponding to a mass rate 
of 10"^ - 10"^ Mq yr-^ which is significantly below the 
observed accretion rates of a few 10~^ Mq yr~^. Hence, we 
suppose that the dynamical influence of the density thresh- 
old should be negligible and its application is acceptable in 
order to avoid very small timesteps and in consequence sig- 
nificantly higher computational costs. 

It is not clear to what extent our res ults are affected 
by the chosen density profile p{r) oc r~^"^. Girichidis et al. 
( [2011 1) find that their results strongly depend on the initial 
profile which they attribute to the varying importance of 
turbulence compared to gravity. As we do not include tur- 
bulent motions in our simulations, our results might depend 
somewhat less on the density profile, in particular with re- 
gard to the accretion rates which also vary rather moderately 



in the work of Girichidis et al. (2011) 



4.2 Disc evolution 

As demonstrated in Fig. [7| a transition from early-type, 
large-scale Keplerian discs to sub-Keplerian discs occurs 
around a normalised initial mass-to-flux ratios /i of ^ 10 
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t = 4000yr log(N [cnT^]) t = 4000 yr log(N [crrT^]) 




r/AU r/AU 

Figure 15. Column density of the disc in the runs 26-20, 10-20, 5.2-20 and 26-0.4 (from top left to bottom right) after 4000 yr. Below 
each slice the Toomre parameter Q (solid lines) and magnetic Toomre parameter Qm (dashed lines) for t = yr, 2000 yr and 4000 yr 
plotted against the radius are shown. Within a radius of approximately 500 AU Q is nearly everywhere below 1 for all runs considered 
thus indicating gravitational unstable discs. Except for run 26-20 fragmentation is suppressed due to the influence of the magnetic fleld 
indicated by Qm — 1. In run 26-20 where even Qm is below 1 fragmentation occurs. 
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independent of the initial amount of rotational energy. This 
is in good accordance with a series of papers studying the 
evolution of low-mass magnetised disc. While 'He nnebelle fc| 
Fromang|(l2008 ) and Hennebelle & Ciardi (2009) find a value 



for fi between 5-10, below which Keplerian disc formation is 
suppressed, [Allen et al.] ( |20Q3| ) and [Mellon Li] (|2008|) find 
no Keplerian discs for up to at least 10. [Puffin Pudritz] 
(2009) studying the possible fragmentation of magnetised 
discs with an initial /i = 3.5 find sub-Keplerian rotation pro- 
files as well in agreement with the results mentioned before. 
Although these simulations apply to low-mass star forma- 
tion with core masses around 1 Mq , the observed maximum 
value of fi, for which the formation of large-scale Keplerian 
discs is suppressed, agrees remarkable well with our finding 
of 5 < /i < 10. 

For the sub-Keplerian discs observed in our simulations 
the infall velocities are generally significantly larger than 
the rotation speed, i.e. Vrot/'^rad < 1, and close to free-fall 
(see right panel of Fig. [4] and Fig. [7|. This is attributed 
to the highly gravitationally unstable configuration of the 
cores containing about 56 Jeans masses. Hence, they cannot 
be stabilised against gravitational collapse by thermal pres- 
sure alone in contrast to low-mass cores containing only ^ 
1 Jeans mass. Consistently, for highly magnetised low-mass 
cores [Hennebelle &: Fromangj ( |2008 ) and [Puffin &: Pudritz 
{"2009 ) find smaller infall velocities as well as t'rot/^'rad ^ 1 
contrary to our results. 

For /i > 10, we find large-scale Keplerian discs in our 
simulations. This is probably a consequence of not having 
turbulence in our initial setup. Turbulence would most likely 
hamper the early formation of large rotating structures and 
delay it to later stages. However, in recent years there is an 
increasing number of observations of rotationally supported 



discs in massive star forming regions (e.g. Cesaroni et al 



2007 and references therein). As an example, discs with 



sizes of a few 100 AU and masses between 0.15 



Fuller 



et al.[[2001[ ) and about 10 Mq ( [Shepherd et al.[[200l| have 
been observed. These discs are similar in size to the discs 
obtained in our runs with weak magnetic fields and also 
show Keplerian-like rotation profiles. In addition, our disc 
masses of 1 Mq calculated in run 26-20 and run 26-4 (see 



Fig. 14) fit perfectly in the mass interval spanned by these 



observations. However, the protostars observed have masses 
around 5-10 Mq and are therefore a factor of a few more 
massive than our most massive sink particles. We attributed 
this to the fact that the observed disc/star-systems are in 
much later evolutionary stages than our systems. 

In contrast, for /i < 10 sub-Keplerian discs are observed 
in our simulations. Again, we mention that our simulations 
end 4000 yr after sink particle formation, thus in a very 
early phase. An observer looking at such a system from edge- 
on will observe a fiattened structure similar to a Keplerian 
disc but without the typical signatures of rotation (see right 
panel of Fig. |9|. Indeed, there is a growing number of ob- 
servations of such fiattened structures reported in literature 
Cesaroni et al.[[2007 and references therein), although 



these observations often refer to more massive structures and 
more evolved protostellar objects than present in our sim- 
ulations. Pue to the insufficient centrifugal support, these 
structures are not in equilibrium and show considerable ra- 
dial infall motions, in accordance with our findings (compare 
bottom panel of Fig. |7|. 



The reason for the lack of centrifugal support in such 
sub-Keplerian discs is the strong magnetic torque acting on 



the midplane (see Fig. 11). Angular momentum is removed 
by magnetic braking at roughly the same rate as it is added 
by the infalling gas. The angular momentum removed from 
the midplane is partly deposited in the outfiow and partly in 
regions further out which are connected to the inner parts by 
the magnetic lever arm ( [Allen et al.| |2003) created through 
the equatorial pinching of the magnetic field (see Fig. [9]). A 
possible way trying to reduce the magnetic braking efficiency 
would be to include non-ideal MHP effects like ohmic dis- 
sipation or ambipolar diffusion in the simulations. However, 
recent numerical work including ambipolar diffusion ( [Mellon 
&: Li|2009[ [Puffin k Pudritz^2009J , ohmic dissipat ion ( |Papp 
&: Basu|2010 ) and also both effects ( Li et al.|2011 ) show that 
even in the case of non-ideal MHP it is not possible to form 
Keplerian discs in such early stages. In fact, these authors 
find that the formation of rotationally supported discs in the 
case of strong magnetic fields is suppressed down to scales 
well below our resolution limit of roughly 10 AU. Hence, at 
the evolutionary stage considered in this work, we cannot 
expect to resolve a proper Keplerian disc even by including 
the effects of ambipolar diffusion or ohmic dissipation in our 
calculations. 

The question now arises how Keplerian discs on 100 AU 
scale, frequently observed around massive protostars, can 
form if the cores have mass-to-ffux ratios ^ 5. The situa- 
tion even intensifies as observed star forming regions usually 
have mass-to-flux ratios which are only slightly supercritical 
with a mean 5 (e.g. Falgarone et al. 2008; Girart et al. 
2009; Beuther et al. 2010 ). Following the long term evolution 



of highly magnetised low-mass cores, Machida et al. (2011) 
find large-scale Keplerian discs occurring after roughly 10 
yr. The authors argue that in their case magnetic braking 
only redistributes the angular momentum within the col- 
lapsing cores but does not remove it. Additionally, they find 
that the outfiows are too weak to ultimately remove a sig- 
nificant fraction of mass and angular momentum from the 
cores. Thus, most of the cores mass and angular momen- 
tum finally fall onto the midplane resulting in large-scale 
Keplerian discs even for high magnetic field strengths. How- 
ever, we note that the authors might underestimate the total 
amount of mass and angular momentum transfered into the 
interstellar medium as only low velocity outfiows are mod- 
eled in their simulations. Considering high velocity jets and 
radiation-driven outflows, the amount of material ejected 
into the interstellar medium could be considerably larger, 
thus impeding the formation of large Keplerian discs. For 
massive stars we expect this to be even more severe as here 
radiation-driven outflows are even more powerful and will 
be able to eject angular momentum - deposited in the enve- 
lope by magnetic effects - at even higher rates than low-mass 
protostellar outffows (e.g. [Arce et al.|2007| . 

At the same time outffows also provide a possible so- 
lution of the disc formation problem by diluting the enve- 
lope in which the magnetic field lines are anchored (see also 
§6.2.2 in [Mellon Li[|2008[) . As t he magnetic braking time 
( [Mouschovias k Paleologou|1980 ) 



Pd 



^A,env Penv 



-1/2 



(13) 



increases with decreasing envelope density, this allows large- 
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scale, centrifugally supported discs to form in particular 
at later stages not considered in this work when most of 
the envelope has been blown away by powerful outflows. 
This would agree with the fact that most of the observed 



However, our observed accretion rates agree well with 
accretion rates from a number of massive star formation 
simulations. To begin with, our accretion rates match those 
necessary to overcome radiation pressure deduced from 1- 



disc/star-systems (see Cesaroni et al. 2007[ and references 
therein) are in a later evolutionary stage than our systems al- 
though the observations could also be biased by the fact that 
massive Class objects are more difficult to detect. The pic- 
ture of a successive growth of centrifugally supported discs 
during the evolution into Class I / II is also supported from 
the theoretical side by Dapp & Basu (2010). Therefore, we 



dimensional calculations (Kahn 19741 Wolfire & Cassinelli 



1987 ). Radiation- hydrodynamical collapse simulations in 2 
dimensions ( |Yorke k Sonnhalter| | 2002| [Kuiper et al.||2010 ) 



and 3 dimensions (Krumholz et al. 2007, 2009; Kuiper e t al.| 

few 

M q yr ~' 1^"^ 



2011 ) with similar core masses reveal accretion rates of c 

A/r„ ,.^-1 io~^ Mq yr~^ very similar to ours. 



expect large-scale Keplerian discs to form in our runs as well 
if the evolution would be followed over a much longer time, 
thus relaxing the problem of catastrophic magnetic braking. 

In summary, our simulations suggest that Keplerian 
discs do not form around massive protostars in the very 
early stage except for unusually weak magnetic fields. Nev- 
ertheless, we expect that centrifugally supported discs will 
build up during the subsequent evolution of the collapsing 
cores. 

Another question not addressed so far is the problem 
of massive binary formation. In neither of our simulations 
do we see indications for massive binaries so far which of 
course could be a consequence of the early stage the sim- 
ulations end. Nevertheless, our results indicate that mas- 
sive binaries possibly form in later evolutionary stages and 
that the initial mass ratio should be far from unity with the 
more massive star sitting in the centre. Observations, how- 
ever, reveal that a significant fraction of massive binaries 



10"^ Mq yr~' up to 1O_M0 yr" 

Peters et al. ( 2010a|b 2011) simulating the long term evo- 



lution of H ll-regions around massive protostars find similar 
accretion rates as well. This suggests that in our simula- 
tions accretion would continue even if the effect of radia- 



tion would be included. Girichidis et al. (2011 ) studying the 
effect of different initial conditions on massive star forma- 
tion find accretion rates around 10~^ Mq yr~^, somewhat 
higher than ours possibly caused by the omission of rotation 
or magnetic fields counteracting gravity and accretion. Sim- 
ilar accretion rates were also found by |Banerjee Pudritz| 
( |2007 ) studying the very early evolution of a collapsed cloud 
core with initial magnetisation and rotation using a signifi- 
cantly higher resolution than in our work. |Hennebelle et alT] 



(|2011[) observe accretion rates of the order of 10 ^ - 10 ^ 



consists of components with nearly equal masses (e.g Mason 



Mq yr~ somewhat smaller than ours which we attribute 
to the larger core size of 1 pc used in their setup. 

Our accretion rates also agree quite well with theoret- 
ical estimates and observations. Calculating the accretion 
rates with the formula given in the theoretical work of |Mc-| 
|Kee Tan ( 2003| ) adapted to our setup we find a value of 



et al.||1998| [Pinsonneault Stanek||2006t . T his apparent about 3 • lO""" Mq yr ^ very similar to the actually observed 



contradiction could be resolved by the work of |Artymowicz| 
k Lubow|(|1996| andlBate k Bonnelll (|1997|. These authors 



simulating circumbinary discs propose that the masses of 
binary components, even when unequal in the beginning, 
tend to equal as accretion within the disc occurs preferen- 
tially onto the lower-mass companion. Initial density pertur- 
bations not consider here could also be a natural way to from 
massive binaries. However, there are indications that frag- 
mentation induced by initial density perturbations might be 
suppressed by magnetic fields (|Hennebelle k Teyssier|2008 ) . 
Thus, different binary formation scenarios might have to be 
at work. 



4.3 Accretion rates 

As shown in Fig. |12| the accretion rates of the different 
runs vary only by a factor of 3 which is attributed to the 
two competing effects of the magnetic field namely magnetic 
braking and the Lorentz force. Adopting accretion rates of 
a few lO"'^ Mq yr ^ as observed in our simulations, stars of 
about 30 Mq would form within some 10^ up to a few 10^ yr 
roughly independent of the initial conditions. However, this 
only holds if the accretion rates stay roughly constant over 
the whole formation process which might be an oversimpli- 
fication as indicated by the slight decrease seen in Fig. [13] 
(see also|Klessen"2001':'Schmeja k Klessen| l2Q04 ) . A possi- 
ble way to significantly change accretion rates would be to 
vary the initial density profile and mass of the molecular 
cloud core, parameters which are not explored in this work 



accretion rates. Observational results for accretion rates in 
massive star forming regions are hard to obtain and are often 
calculated only indirectly via observed outflow mass rates. 
Nevertheless, results from several high-mass star forming re- 
gions indicate accretion rates of the order of 10~^ - 10~ ^ Mq 
yr"^ (e.g. Beu ther et al.|2002a|c| [Beltran et al.||2006| ) again 
in good accordance with our results. 

As shown in Fig. [13] the accretion rates do not drop 
significantly over time (except for runs with /x = 2.6, see 
Section 4.1 ). Hence, the outflows launched seem to be not ca- 



in order to limit the computational costs (see e.g. Girichidis 
eraL][20TT| . 



pable of significantly reducing mass accretion onto the sinks 
over time. In order to analyse to what extent the combined 
effect of magnetic fields and rotation influences mass accre- 
tion, we performed the reference calculation inf-0 with no 
magnetic field and zero rotation. In this run the sink parti- 
cle accretes 3.55 Mq within the first 4000 yr corresponding 
to a time averaged accretion rate of 8.86 • 10~^ Mq yr~"^. 
Hence, the accretion rates in runs with magnetic fields and 
rotation are reduced to a level of about 35% - 95% of this 
value (see Table [2|). Our simulations show that feedback due 
to outflows is not able to stop accretion onto the protostars 
and hence outflows are not able to determine the low star 
formation efficiency observed in molecular clouds. 



5 CONCLUSION 

We have studied the collapse of massive molecular cloud 
cores with varying initial rotational and magnetic energies. 
The cores are supercritical with mass-to-fiux ratios between 
2.6 and 26 and have rotational energies well below the gravi- 
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tational energy. Containing about 56 Jeans masses the cores 
are highly gravitationally unstable and hence might be sites 
of protostellar cluster formation. We focussed our discus- 
sion on the formation of protostellar discs and on protostel- 
lar accretion as measured by sink particles. We find that 
disc properties are highly sensitive to the initial magnetic 
field strength whereas protostellar accretion rates are only 
marginally influenced by varying initial conditions. In the 
following we summarise our main findings. 

1. For normalised mass-to- flux ratios fi below 10 the for- 
mation of centrifugally supported discs is completely sup- 
pressed. Instead, for /x < 10 sub-Keplerian discs are formed 
showing nearly no rotational support but radial infall veloc- 
ities close to free-fall. For weak magnetic fields (/i > 10), 
however, well defined Kepler ian discs with sizes of a few 100 
AU build up over time. This agrees well with several studies 
of collapsing low-mass cores. 

2. Sub-Keplerian rotation for strong magnetic fields (/i < 
10) is caused by magnetic braking. Analysing the torques 
acting on the midplane we find that angular momentum is 
removed due to magnetic braking at roughly the same rate as 
it is added due to the gas infall thereby preventing Keplerian 
discs from forming. 

3. Observed accretion rates are of the order of a few 10~^ 
Mq yr~^ varying only within a factor of ^ 3 between the 
individual runs. This variation is remarkably small consid- 
ering the large differences in the initial conditions varying 
over more than two orders of magnitude. We attributed this 
to two competing effects of the magnetic field. Increasing 
the magnetic field strength results in an increased accre- 
tion rate due to an enhanced magnetic braking efficiency 
lowering centrifugal support. Above a certain field strength, 
however, a further increase leads to declining accretion rates 
due to the effect of magnetic pressure and tension. This re- 
sults in rather moderate changes in the accretion rates for 
different initial conditions. Furthermore, accretion rates for 
different amounts of angular momentum converge with in- 
creasing field strength due to the effect of magnetic braking. 

4. For the majority of the simulations disc fragmentation 
does not occur. Analysing the Toomre parameter Q in these 
discs reveals that fragmentation is mainly suppressed due to 
the magnetic pressure. Thermal pressure alone cannot ac- 
count for the stabilisation as Q is well below 1 whereas the 
magnetic Toomre parameter Qm settles around 1 indicating 
stability. In two simulations, which are subject to fragmen- 
tation, neither thermal nor magnetic pressure can stabilise 
the discs. Accordingly, both Q and Qm are well below 1. 

5. In the two runs with disc fragmentation more than 10 
sink particles are formed during the first 4000 yr. Among 
these sinks only the first one created reaches a mass above 1 
Mq thus containing more than 80% of the totally accreted 
mass. All other particles have masses well below 0.1 Mq. 
The discs formed in both cases reach masses around 1 Mq 
somewhat below the totally accreted sink particle masses. 

6. The outflows launched from the protostellar discs are 
not capable of significantly reducing mass accretion over 
time. Compared to a run with zero magnetic field and zero 
rotation, mass accretion is reduced to a level of 35% - 95%. 

7. Radial profiles of column density and temperature ex- 
hibit accretion shock features moving outwards as time 
evolves. The shocks occur when the infalling gas hits either 



centrifugal or magnetic barriers thus strongly decreasing the 
infall speed. 

A growing number of observations of discs and bipolar 



outflows around high-mass protostars (see Beuther & Shep- 
herd 2005 Cesaroni et al.|l2007, for recent reviews) sup- 



port a high-mass star formation scenario via disc accretion. 
On the other hand, observations also reveal that prestellar 
cores with masses ranging from 2 - 2000 Mq are usually 
only slightly supercritical with < 5 ( Falgarone et al.|2008 



Girart et al.||2009; Beuther et al.||2010). Together with our 



numerical results this suggest that there should be no Kep- 
lerian discs in the very early stages of typical high-mass star 
forming regions but rather flattened, strongly sub-Keplerian 
structures. This apparent dichotomy has an important im- 
pact on the formation of discs around massive stars. To en- 
able the observed presence of centrifugally supported discs 
in later stages, effects in the later evolution of the system 
are required reducing the efficiency of magnetic braking. We 
discussed possible effects in the context of massive star for- 
mation like outflows and non-ideal MHD leading to a suc- 
cessive growth of discs in the later evolution. 

So far, the influence of turbulence is completely ne- 
glected in our setup. Hence, we plan to include initial veloc- 
ity perturbations in our setup to obtain more realistic initial 
conditions as indicated by observations (e.g. Caselli & Myers 
1995| ) . For subsequent simulation it would be also interesting 
to study the effect of ambipolar diffusion which is expected 
to act efficiently in the high density regime of protostellar 
discs. For this purpose, the existing set of simulations serves 
as a useful guide to select representative simulations to be 
repeated with increased resolution and additional physics. 
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Table 3. Performed simulations for the resolution study with 
maximum spatial resolution, threshold density for sink particle 
creation, formation time of the first sink particle and total ac- 
creted mass after 2000 yr. 
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(kyr) 


(Mq) 


26-4-L 


18.9 
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15.1 
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26-4 
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15.3 
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APPENDIX 

Here we compare simulation 26-4 to two more runs with 
identical initial conditions but a maximum spatial resolu- 
tion varied by a factor 4 in either direction. In particular, 
the initial resolution is identical to that in run 26-4, i.e. in 
the beginning the mesh in the core has a spacing of 302 
AU. We list the runs and their corresponding parameters 
in Table [3] The critical value of the density above which 
sink particles are created is adapted in accordance with the 
resolution. For all runs performed the refinement criterion 
applied guaranties that the midplane region is resolved on 
the highest level used. In particular, we focus in this analysis 
on accretion properties and radial profiles of different quan- 
tities in the midplane. Due to computational cost reasons 
run 26-4-H is followed for 2000 yr only. Hence, we compare 
the results at this time. 

First, we consider radial profiles of the density, tempera- 
ture and velocity in the midplane. The quantities in each run 
are taken in the midplane. In the density profile (left panel of 
Fig. 16) the accretion shock occurring at ^ 150 AU is clearly 
resolved in the runs 26-4-H and 26-4. In run 26-4-L, how- 



ever, the shock is somewhat smoothed out due to the limited 
resolution. Hence, a resolution of 4.7 AU seems required to 
properly resolve the accretion shock. Within the accretion 
shock, however, the density increases with resolution. We 
attribute this to the fact that the vertical structure of the 
disc is not fully resolved at least in the runs 26-4 and 26-4-L. 
Here in large parts the vertical disc height is represented by 
a few grid points only (compare Fig. [9]). Therefore, to fully 
resolve the vertical disc structure, a higher resolution, prob- 
ably even above that in run 26-4-H, would be needed which 
is not feasible for computational costs reasons. 

A similar result holds for the temperature profiles as 



well (middle panel of Fig. 16 ). The temperature jump at the 
accretion shock seems be to reasonably well resolved in run 
26-4 whereas in run 26-4-L it is smoothed out markedly. 
Within the accretion shock the temperature reaches as 
higher values as higher the final resolution is. This is due 
to the strong coupling of temperature and density above 
10~^^ g cm~^ where the gas gets optically thick resulting in 
higher temperatures at higher gas densities. 

Next, we analyse the velocity structure in the midplane 




(see right panel of Fig. 16). In run 26-4-L no Keplerian disc 
has built up yet. However, we mention that in its further 



2000 



Figure 17. Total accretion rate for the first 2000 yr in the runs 
26-4-L (blue hue), 26-4 (black line) and 26-4-H (red hue) with 
a maximum spatial resolution of dx = 18.9 AU, 4.7 AU and 1.2 
AU, respectively. The mean accretion rates (see also Table |3j are 
as higher as lower the resolution is but seem to converge towards 
higher resolution. 



evolution the rotation reaches Keplerian velocities as well. 
In contrast, in run 26-4 the rotation is already Keplerian up 
to a radius of ^ 60 AU in good agreement with run 26-4- 
H. Furthermore, the comparison shows that in run 26-4 the 
decline in Vrot/'^^kep at 20 AU is most likely a resolution ef- 
fect as in run 26-4-H this decline occurs at a roughly three 
times smaller radius of ^ 8 AU. This supports the state- 
ment made in Section [3. II that the inner 10 AU in runs with 
a resolution of 4.7 AU are strongly affected by numerical res- 
olution (see also Fig.|7|. The radial velocities at radii larger 
10 AU show a reasonably well agreement between all runs. 
Hence, we conclude that regarding the velocity structure in 
the midplane run 26-4 is reasonably well converged. 

In Fig. ^1 we show the time evolution of the total ac- 
cretion rate for the three runs considered. In run 26-4-H two 
more sink particles are created after roughly 1500 yr which 
cause the large variations in the accretion rate. In general, 
however, the accretion rates of the three runs are of the 
same order of magnitude. From the formation time to listed 
m Table [3] it can be inferred that to increases with spatial 
resolution. This behaviour is expected as for higher resolu- 
tion sink particles are created at higher densities and thus 
later times during the collapse. It can also be inferred from 
Table |3] that the mean accretion rate, i.e. the total accreted 
mass, decreases with increasing spatial resolution. The ac- 
cretion rate of run 26-4-L is higher than that of run 26-4 by 
about 35%. We therefore conclude that regarding accretion 
properties run 26-4-L is not yet converged. The difference 
in accreted mass between run 26-4-H and run 26-4 is of the 
order of 2% thus significantly lower than the difference be- 
tween run 26-4-L and run 26-4. Hence, there seems to be 
a convergence of the accretion rates with increasing resolu- 
tion and, although run 26-4 seems not fully converged, we 
conclude that a resolution of 4.7 AU is sufficiently high to 
properly describe accretion properties of the protostars. 

In summary, we can say that a resolution of 4.7 AU 
seems sufficiently high to properly follow protostellar accre- 
tion rates, resolve the accretion shock at the edge of the disc, 
and display the velocity structure in the disc down to about 
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Figure 16. Radial profile of density (left), temperature (middle) and velocity structure (right) after 2000 yr for the runs 26-4-L (blue), 
26-4 (black) and 26-4-H (red). The accretion shock seen in the left and middle panel is not well resolved for run 26-4-L. 



10 - 15 AU. In contrast, run 26-4-L with a resolution of 18.9 
AU reveals significant differences from run 26-4-H showing 
that is not yet converged. However, also run 26-4 seems to 
be not fully converged regarding to the density and temper- 
ature structure in the disc. We attribute this partly to the 
poor resolution of the vertical disc structure. Hence, a higher 
resolution, probably even above the one used in run 26-4- 
H, would be necessary to reach convergence regarding this 
point. However, due to significantly higher computational 
costs this has not been feasible wherefore a spatial resolu- 
tion of 4.7 AU is used throughout the paper. This particular 
choice is also motivated physically as we want to resolve the 
first core. Hence, the Jeans length at densities above 10~^^ 
g cm~^ has to be resolved requiring a resolution of a few 
AU. 



